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This paper originates from the joint efforts of an aeroelastic study team in the Applied Vehicle Technol- 
ogy Panel from NATO Science and Technology Organization, with the Task Group number AVT-191, titled 
“Application of Sensitivity Analysis and Uncertainty Quantification to Military Vehicle Design.” We present 
aeroelastic uncertainty quantification studies using the SemiSpan Supersonic Transport wind tunnel model at 
the NASA Langley Research Center. The aeroelastic study team decided treat both structural and aerodynamic 
input parameters as uncertain and represent them as samples drawn from statistical distributions, propagat- 
ing them through aeroelastic analysis frameworks. Uncertainty quantification processes require many function 
evaluations to asses the impact of variations in numerous parameters on the vehicle characteristics, rapidly in- 
creasing the computational time requirement relative to that required to assess a system deterministically. The 
increased computational time is particularly prohibitive if high-fidelity analyses are employed. As a remedy, 
the Istanbul Technical University team employed an Euler solver in an aeroelastic analysis framework, and 
implemented reduced order modeling with Polynomial Chaos Expansion and Proper Orthogonal Decomposi- 
tion to perform the uncertainty propagation. The NASA team chose to reduce the prohibitive computational 
time by employing linear solution processes. The NASA team also focused on determining input sample distri- 
butions. 


Nomenclature 


LHS Latin Hypercube Sampling 

MCS Monte Carlo Sampling 

PCE Polynomial Chaos Expansion 
POD _ Proper Orthogonal Decomposition 
ROM _ Reduced order model 


I. Introduction 


eroelasticity, as a multidisciplinary research field, investigates the behavior of an elastic structure in the airstream 
Aaa the interactions between inertial, aerodynamic, and structural forces. With the use of advanced technologies 
in aircraft design and analysis, the need for predicting such aerostructural interactions by using high-fidelity compu- 
tational models has increased. Most analysis procedures involving high-fidelity modeling consider the system and 
solution as deterministic. However, as computational models have advanced, parametric sensitivities have become 
more intrinsic and more hidden. Recent work has demonstrated solution sensitivity to small changes in geometry, 
freestream conditions and computational parameter specifications. This sensitivity leads to considering these simu- 
lations from an uncertainty standpoint. As another motivation in experiments, input parameters are only controlled 
to within some precision level and the input and output parameters are only known to within measurement errors. 
Comparisons between computations and experiments should consider these uncertainties. 
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In the current paper, we endeavor to recast our aeroelastic simulation inputs and outputs to treat the analyses as 
stochastic processes. Two different approaches and philosophical viewpoints are utilized to quantify the effects of 
uncertainties in aeroelastic systems. The methods, modeling and results are produced from research efforts of teams 
from NASA and Istanbul Technical University (ITU). The NASA work utilizes linear aerodynamic modeling, while the 
ITU work utilizes higher fidelity aerodynamic modeling and a reduced order model (ROM) approach for propagating 
the uncertain parameters. 

A flexible vehicle configuration, the NASA Semispan Supersonic Transport (S*T) wind tunnel model, was chosen 
as the basis for uncertainty studies and propagation. The current work utilizes information, analytical models and data 
reported by the (S*T) project team in 2012 at a special session at the 53rd AIAA/ASME/ASCE/AHS/ASC Structures, 
Structural Dynamics and Materials Conference. A program overview paper was authored by Silva*® that documents 
the wind tunnel model, four wind tunnel tests conducted between 2007 and 2010, and references laboratory tests and 
analyses. The products from those efforts are utilized as the starting point for the current uncertainty analyses; these 
aspects will be further identified in the body of this paper. 

The current paper focuses on defining statistical distributions of input parameters, propagating them through the 
two analysis processes, and examining the statistical distributions of output parameters. The analysis approaches will 
be described first, followed by a description of the development of an example input parameter distribution to represent 
uncertainty associated with a structural parameter. Example results from each of the analysis methods will then be 
presented, followed by discussion and concluding remarks. 


II. Analysis Approaches 


A. Uncertainty Propagation Using Linear Analyses 


The methods, modeling and results generated by the NASA team were produced from linear aerodynamic and struc- 
tural computations. The models used as the baseline cases came from the project team that designed, tested and 
analyzed the configuration. The major steps taken in performing the uncertainty propagations were: 

e Define input variation 


e Generate random input samples 

e Generate input file for each input sample 

e Perform simulation for each input sample 

e Extract output parameters from each simulation 


e Perform statistical analyses of the evaluation parameters 


Three types of simulations were performed in the NASA analyses: static aeroelastic analyses, normal modes anal- 
yses and flutter solutions. Each of the analyses were performed using standard analysis options in MSC NASTRAN.!® 
Input and output parameter screening was performed using static aeroelastic analysis. Normal modes analyses were 
performed to show the sensitivity of modal frequencies to the structural variations. Flutter analyses were performed to 
show the sensitivity to modal parameters in addition to sensitivity to the more fundamental input parameters. 


B. Reduced Order Modeling Using Euler Simulations 


The ITU team performed static aeroelastic and flutter analyses of the S*T wing by using the Polynomial Chaos Expan- 
sion (PCE)??~*> and Proper Orthogonal Decomposition (POD)?”-29-3!.33:34 based Reduced Order Models,”°3? which 
were generated by using an Euler flow based aeroelastic solver in ZEUS! software. The aeroelastic solver of ZEUS 
employs mode shapes obtained from the NASTRAN modal solver. The ROM results were compared to the initial com- 
putational results of the high-fidelity models. These ROMs could predict the responses in a small input data range, and 
in this case both the PCE and POD methods gave satisfactory results while their accuracies were thoroughly discussed 
in the present work. Secondly, both the PCE and POD methods were used to construct a ROM, which was accurate 
when large variations of input parameters were considered. The flutter behavior of the ST wing could accurately be 
predicted over the entire flight regime for Mach numbers ranging from 0.60 to 1.20 by the PCE and POD based ROMs. 
This will be demonstrated in this paper through detailed comparisons of the flutter results of the ROMs to the initial 
high-fidelity computational analyses. 

Uncertainty distributions of structural and aerodynamic input parameters were then propagated using the ROMs. 
Output parameters from the reduced order static aeroelastic and flutter analyses were examined to assess the influences 
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of the input parameters that were varied. The input sampling was performed using Monte Carlo Simulation (MCS) and 
Latin Hypercube Sampling (LHS) techniques, which are frequently used to propagate uncertainties in a probabilistic 
framework. The input parameters that were treated as uncertain in these analyses were angle of attack, Young’s 
modulus and dynamic pressure. The influence of the first two were assessed for both the static aeroelastic and flutter 
analyses. The dynamic pressure uncertainty was only evaluated for the static aeroelastic cases. The output parameters 
that were examined were flutter frequency, flutter speed, flutter dynamic pressure, lift coefficient (Cz), drag coefficient 
(Cp), pitching moment coefficient (C,,), and wing tip twist angle. 


I. Polynomial Chaos Expansion (PCE) Methodology 


Nonintrusive PCE defines a ROM basis in terms of orthogonal polynomials. In the present study, Hermite polynomials 
are used as the orthogonal polynomials since the input parameters are assumed to vary with respect to Gaussian 
distribution. The definition of the first few Hermite polynomials are given below in terms of the standard variable, &: 


H,(§) =1 (1) 
M(E)=§ (2) 
(§)=6?-1 (3) 
H3(§) = 6° —36 (4) 
H4(§) = 6° - 667 +3 (5) 
Hs(§) = &° — 108° + 156 (6) 
_ £6 4 2 
H6(§) = 6° — 156" + 45€- — 15 (7) 
The chaos coefficients can be calculated by using the following relation for a general physical system: 
uy Ajo Ay Ay. . Ap ay 
uz Ayo Hx, Hx. . Hap a2 
- (8) 
Un Ano An An. . Anp ap+1 


where u; (i = | to n) is an output response that is known from the analyses and a; (j = 1 to p+ 1) shows chaos 
coefficients to be determined. The order of the approximation is denoted by p and n shows the dimension. In this 
matrix system, the parameters H;; indicate the multiplication of H; and H; functions defined in Eq. (4) to (10). 
However, in the two-dimensional case, H; is the function of the first standard random variable, &, while H; is the 
function of the second standard random variable, 7. Since three input parameters are defined in the present work, 
Hijjx = Hi(§)-Hj(n) - Hy) (where y is also a standard random variable) multiplications should be computed. The 
ROM can be constructed by using the chaos coefficients. In order to calculate the chaos coefficients, linear regression 
can be applied. Considering the standard normal variables €, n and y, a general linear regression model, U, can be 
defined in closed form as below: 


U=Ha+e (9) 


where € denotes the error between real values and approximate values. The regression coefficients, @, can be computed 
by using Eq. (10): 


a4=(H'H)'H'U (10) 


After determining the regression coefficients, the ROM to represent the reference initial design can be constructed for 
the two input variable case by using the relation below: 


Uj(§*,n") = a Hjo(S*,n") +a2H i (S*,1*) +... +ap41Ajp(S",n") (11) 


where &* and 1* are the values of the standard random variables for the case that needs to be computed. 
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2. Proper Orthogonal Decomposition (POD) Approach 


The Proper Orthogonal Decomposition (POD) method is used to generate a ROM. The system is defined as a linear 
combination of independent orthogonal basis functions. The orthogonal functions, called the “POD basis functions,” 
are formulated such that the effect of the first mode is the highest and the second mode is less effective and so on. A 
physical model with n outputs can be expressed as a linear combination given as below: 


Dy] (nx1) = aj{® {| (nx1) (12) 


Ms 


1 


Pi 


where y is an output response vector for a snapshot of the input variables, a; are POD coefficients, ®; are POD basis 
vectors and k is the order of ROM (k <n). In this study, the “method of snapshots”*> is used for creating ROMs. 
A snapshot matrix T, is an n X m matrix (m is number of snapshots) and can be constructed by using the output 
parameters of the system: 


IT] (nxm) = 


yh) : . yy” 


(nxm) 


Any physical system given with a snapshot matrix can be represented by a surrogate response function which super- 
poses the basic modes of the model. The correlation matrix is defined in terms of the snapshot matrix and the number 
of POD basis vectors: 


1 
[CV (uxm) = = (ux) * TVmxn) (13) 


The POD basis vectors and coefficients are defined through the eigenvalues and eigenvectors of the correlation ma- 
trix [C]}® = A®. Singular Value Decomposition (SVD) facilitates the construction of the POD based ROMs. In this 
method, the factorization of a real matrix, T is given below: 


LE lesen = [)nxn) FE lvsenn) VV" xm) (14) 


where U, the columns of which consist of the eigenvectors of [T] (T}", is the left singular vector of T and these 
eigenvectors are POD basis vectors. V, the columns of which consist of the eigenvectors of (T)" [T], is the right singular 
vector of T. [X] is a diagonal matrix representing the square roots of POD eigenvalues, A;. The POD coefficients are 
calculated by projecting the snapshots on the POD basis. 


aj= TH” (1x) © [Pylinx) (15) 
The number of POD basis vectors is equal to the number of snapshots m (and also POD coefficients), however, the 


number of retained POD basis vectors (POD modes) k is determined by considering the captured energy from the 
original model. The total captured energy is defined as: 


k 
Etotal = >. Xi (16) 
i=1 


where A is the vector of normalized eigenvalues of the correlation matrix, C, and k is the required number of POD 
modes to reach the specified captured energy level. In the present work, the captured energy level is taken as 99.9 % 
in all applications. 
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III. Development of Statistical Distributions of Input Parameters 


Development of statistical distributions to represent input parameter variations were performed by the NASA 
team. The information required to formulate the distributions of structural parameters is discussed in this section. 
The model hardware, testing and finite element model (FEM) tuning process were the key elements in developing the 
input parameter distributions. The variations for other input variables that are utilized in the current work—angle of 
attack, dynamic pressure, Mach number and modal frequencies—are not discussed in detail here; the processes and 
assumptions for developing these additional input parameter variations are detailed in the NATO-STO-AVT-191 Task 
Group Report.!7 


A. Hardware Characterization Testing 


The S*T wing is a complex aeroelastic semispan wind tunnel model, which is designed for aeroelastic/ aeroservoelastic 
analysis.*” The fuselage consists of a graphite-epoxy flexible beam attached to an aluminum C-channel rigid beam. 
Wingbox construction consists of an upper and lower skin, each made from fiberglass/epoxy face sheets sandwiching 
a paper honeycomb core. Figure | shows the configuration. 


Active all-moveable 


: : Active ride 
horizontal tail 


Control vane 


Veo: | Flexible wing | 


Engine nacelles w/ 
Flexible attachments 


Figure 1. S*T configuration. 


Extensive wing geometry measurement, stiffness testing and modal survey testing (ground vibration testing, GVT) 
were performed. These results were used to improve the aerodynamic surface definition and the structural model. 
The flexible fuselage beam element material properties were used as principal tuning parameters during model update 
procedures. Figure 2 shows the results from tuning the analytical model to match the measured deflections. The range 
of the modifications used in performing this tuning serve as the basis for the material property distributions used in the 
uncertainty propagation studies. 

The GVT results*’contain measured values for the frequencies and damping values associated with the first six 
modes. The analytical and measured frequencies and descriptions of the mode shapes are shown in Table 1.°° Compar- 
isons of the analytical and measured frequencies give indications of the success level of the stiffness tuning exercise. 
The scatter among the measured values for each modal frequency also provides an indication of the variation associated 
with measuring the frequencies. Measured damping values, reported in reference 37, contained significant scatter. 


B. Structural Model 


The FEM utilized in this study for the baseline analysis cases is the unmodified model from the S*T team that is 
detailed by Sanetrik.** The FEM elements, emphasizing the fuselage, are shown in Figures 3, 4 and 5. Figure 5 in 
particular shows the flexible fuselage beam material property segments that will be modified as part of the uncertainty 
propagation studies. The numbers given in the figure correspond to the designations used for material property charac- 
terization sets. Associated with each material property set is a value for density, Young’s modulus, the shear modulus 
and Poisson ratio. In the current studies, the only material property to be treated as an uncertain variable is the Young’s 
modulus. 

The fuselage beam elements shown in Figure 5 were considered the primary candidates for parametric variation in 
this study, based on the tuning processes of the S*T team. The input variable ranges were based on the values used by 
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Figure 2. Finite element model tuning by S*T test team; Deflection in degrees ys Flow-wise (chordwise) station. 
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Figure 3. Finite element structural model. 


the S*T test team to tune the FEM. During this procedure, the StT team tuned the six beam materials in the baseline 
structural model, correlating with measured deflection data, as shown in Figure 2. In the current study, these same six 
material properties were evaluated, using minimum, nominal and maximum values employed in the original tuning 
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Table 1. Measured and analytical modal frequencies, Hz. 


Modal frequencies, Hz 


Mode Description Analytical nominal value —- Measured values 
1 Pitch about forward nodal mount 6.29 6.395 6.375 6.249 
2 Pitch about balance 8.43 8.089 7.935 7.838 
3 Fuselage Ist bending 10.07 10.323 10.312 10.059 
4 Wing Ist bending 11.4 11.635 - 11.781 
5 - 12.88 12.528 12.585 12.384 
6 - 14.25 16.571 15.066 15.603 
7 - 15.28 - - - 
8 - 17.6 - - - 


+ +444 
++ 
TTI} 


+O 


Forward 


Nodal Mount Wing-to-fuselage Connections Aft Nodal 


Point, Mount Point, 


U-spring Fuselage Beam Elements U-spring 


Figure 4. Finite element model: Connection points from fuselage to wing, nodal mount points and fuselage beam elements. 


Figure 5. Beam elements’ material property numbers. 


process. The values chosen by the S4T team to use in the final version of the FEM were designated as the nominal 
values. 

Three input distributions and drawn sample sets characterizing the uncertainty of the Young’s modulus were gen- 
erated: 1) a normal distribution sampled with LHS, 500 samples drawn without replacement; 2) a normal distribution 
sampled with LHS, 25 regions with 25 resamplings with replacement (625 total samples); and 3) a uniform distribution 
sampled with LHS, 25 regions with 25 resamplings with replacement (625 total samples). For material 21, which will 
be shown subsequently to be of greatest interest, the Gaussian distributions from which samples were drawn had a 
mean value of 6.390e+06 psi and standard deviation of 0.945e+06 psi. The detailed statistics for all three distributions 
for this material property are given in Table 2. 
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IV. Linear Analyses 


A small subset of the linear analysis results generated in this study are discussed and presented in the current paper. 
Additional results and comparisons with different uncertainty quantification approaches are detailed in reference 17. 
In the current section, an example of parameter screening using static aeroelastic analysis is presented, along with an 
example case of propagating a distribution of an input parameter through a normal modes analysis. Finally, example 
results are shown for propagating both uniform and normal distributions of an input parameter through the flutter 
analysis process. 


A. Static Aeroelastic Analyses for Parameter Screening 


Material property screening studies used the flexible-to-rigid ratio of the lift coefficient due to the horizontal tail 
deflection as one of the evaluation parameters. This particular output variable was chosen due to several factors: 1) 
this parameter was shown to be linear with respect to dynamic pressure; 2) the S*T program’s main research focus 
was on aeroservoelasticity, which emphasizes the criticality of control derivatives and 3) most control law designers 
associated with the program chose the horizontal tail as their actuation surface for performing aeroservoelastic tasks. 
The screening analysis results, shown in Figure 6, indicate that the Young’s modulus associated with the material 
property 21 has the greatest influence. This is shown by the relative slopes of each of the lines; the slope associated 
with material property 21 is the steepest. 
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Figure 6. Static aeroelastic analysis parameter screening; Flexible-to-rigid ratio variation of C,,7 due to variation of 
flexible fuselage beam Young’s moduli. 


Thus, the Young’s modulus of material 21 was chosen as the structural model parameter to vary and propagate 
through static and dynamic aeroelastic simulations. This material property is associated with beam elements compris- 
ing the model’s fuselage in the region near the mid-chord wing root, as shown in Figure 5. 


B. Normal Modes Analysis 


The second major type of analysis that was performed using linear methods was normal modes analysis. For this 
analysis, the input parameter identified as the most significant using the screening analysis was varied—the Young’s 
modulus of material property 21. Any variation in dynamic pressure, Mach number, angle of attack or structural 
damping would have no impact on the output parameters of a normal modes analysis. The output parameters that are 
examined here are the frequencies of the first (lowest value) 30 natural modes. The three input distributions shown in 
Table 2 for the Young’s modulus variation were each used as inputs to this process. 

Histograms of the first 30 natural frequencies are shown in the 30 subplots of Figure 7 for the second input case 
as numbered in Table 2. Notice that for the input distribution, which is Gaussian, some of the modal frequency output 
distributions appear Gaussian and some appear skewed (i.e., some of these output distributions are likely better fit by 
beta distribution than by Gaussian distributions). There are also a few frequencies where the influence of the input 
variation is so small that the histogram shows quantization levels associated with the significant digits of the output 
data file (i.e., the first two subplots in the third row, associated with modes 13 and 14). The results for the second 
input case are presented in Figure 8. In Figure 8, the mean values for the first 8 modal frequencies are shown by the 
horizontal lines. The variations in the simulation outputs are captured by the heights of the blue rectangles, indicating 
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Table 2. Statistical analysis of input distribution samples, Young’s modulus varied. 


Case Input Sampling Mean Standard Skewness’ Kurtosis Coefficient 
distribution Method deviation of 
type variation 


1 Gaussian LHS 500 regions, 6.39E+06 9.45E+05 2.39E-02 2.96E+00 1.48E-01 
no replacement 
2 Gaussian LHS 25 regions 6.39E+06 9.49E+05  -5.91E-02 3.04E+00 1.49E-01 
25 resamplings, 
with replacement 
3 Uniform LHS 25 regions 6.39E+06 1.09E+06 -4.20E-03 -2.90E-01 1.71E-01 
25 resamplings, 
with replacement 


Mode #’s 
1-6 


Aa 4 4 
ae ae ee | 
ee allt dh 
ae | a 

y 7 


25 - 30 


Figure 7. Histograms of the first 30 modal frequencies; Results from LHS (25, 25) samples from Gaussian distribution of 
the Young’s modulus of material 21. 


the +/- 3 sigma bounds of the frequencies. The figure shows that only the second and third modes are significantly 
modified by the variations applied to the Young’s modulus. 


C. Flutter Analysis 


Flutter analyses were conducted for eight baseline cases, and for input parameter variations of each of those eight cases. 
The baseline cases correspond to the subsonic Mach numbers where the S*T team generated computational results. 
Three uncertainty distributions were propagated for each condition for five input parameter sets: structural damping, 
Young’s modulus, Mach number and natural frequencies. Each input distribution was propagated independently. The 
complete results are given in reference 17. In the current section of this paper, example results are shown for two 
distributions representing the damping variation. The Gaussian distribution of damping was generated with a mean 
value of 0.02, and a standard deviation of 0.006; the uniform damping distribution was generated with a minimum 
value of 0 and a maximum value of 0.04. 

The doublet lattice method was employed to generate the aerodynamics and the PKNL (p-k nonlinear) flutter 
solution of MSC/NASTRAN was employed. For each of these baseline cases, the structural damping was set to zero, 
and the first thirty flexible modes were included in the structural definition, corresponding to including all modes 
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Figure 8. Variation of Young’s modulus for material 21, Gaussian LHS (25, 25); Influence on the first 8 modal frequencies. 


below 100 Hz. The flutter analyses were performed echoing the baseline process performed by the S*T team, keeping 
the velocity constant and varying the density. Although flutter results are shown here at transonic Mach numbers 
(Mach 0.80 and 0.95), this should not be interpreted as an endorsement to utilize linear methods in this range. The 
results are possibly relevant at Mach 0.8 for low angles of attack, but this must be established by comparisons with the 
experimental data or higher fidelity simulations. 

In linear analysis, the flutter condition is defined as the dynamic pressure that causes a mode to have a damping 
value of zero. The structural damping appears in the flutter solution equations as an additive term. Once the aeroelastic 
modal damping has been calculated, changes due to the structural damping can be approximated by using the structural 
damping as the definition of neutral stability, rather than using the value of 0. This is not an exact calculation due to 
the presence of nonlinearities introduced by the generalized aerodynamic forces, which are calculated as functions of 
reduced frequency and thus are functions of structural damping. These influences are small and are neglected here. 

The derivative of the aeroelastic damping with respect to dynamic pressure is referred to as the steepness of the 
flutter crossing. A steeper flutter crossing corresponds to a case that is more sensitive to aerodynamic damping (or 
aerodynamic energy input) and less sensitive to structural damping. The higher Mach number cases are more sensitive 
to aerodynamic loading than those at lower Mach numbers. Conversely, the lower Mach number cases—with flatter 
flutter crossings—are more sensitive to structural damping than the higher Mach number cases. The result of the 
difference in the flutter crossing steepness at different Mach numbers is illustrated in Figure 9, where the flutter 
dynamic pressure variation is observed to be considerably larger at Mach 0.6 than at the higher Mach numbers. 

We note that for a Gaussian input distribution of the damping, the output flutter dynamic pressure is also Gaussian. 
The damping variation flutter results are plotted with the baseline analysis results (shown by the red circles) in Figure 9. 
In the figure, the output histograms are shown on their sides. The green bar graphs show the histograms of the flutter 
dynamic pressures due to a uniform input variation of the damping. The height of the histogram (in terms of flutter 
dynamic pressure) is greatest for the lowest Mach number indicating that the lower Mach number flutter onset is 
more sensitive to uncertainty in damping than the higher Mach numbers. This reflects the steeper flutter crossing as 
discussed above. The uniform input distributions result in uniform output distributions. The blue bar graphs shows the 
histograms of the flutter dynamic pressures due to Gaussian input variation of the damping. In comparing the results 
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among the Mach numbers, the trend in the range of flutter dynamic pressures are similar to those discussed for the 
uniform variation. The shapes of the output histograms are Gaussian, reflecting the linear relationship between the 
flutter dynamic pressure to the damping. 
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Figure 9. Flutter boundary, including damping variation results. 


V. Euler Based Aeroelastic Analysis and Reduced Order Modeling in UQ 


In this section, an Euler solver is used to compute the static aeroelastic and flutter responses of the S*T wing for 
uncertainty quantification (UQ) process. The Euler solver provides a good balance between the complete modeling of 
the flow physics and the computational speed needed for repeated function evaluations in UQ. The static aeroelastic 
and flutter analyses are accomplished for the given reference cases at different flight conditions with various Mach 
numbers, angles of attack values and dynamic pressures. After performing the sufficient amount of analyses, these 
high-fidelity computational models are used to produce ROMs. These models enable considerable reduction in com- 
putational time while investigating the effects of uncertain input variations on output aeroelastic responses through 
uncertainty quantification. In the present work, the reduced modeling strategy is used to represent the high-fidelity 
computational models of the S*T wing to present accurate solutions while providing a considerable reduction in the 
computational time. 

The fundamental steps of the ROM strategy of this work first require construction and analysis of the high-fidelity 
computational models. These computational models are developed to perform static aeroelastic and flutter analyses of 
the S*T wing in the ZEUS software. ZEUS is an Euler unsteady aerodynamic solver developed for aeroelastic solutions 
of complex geometries.!° It involves an automated mesh generation scheme and overset grid capability for complex 
configurations while using a Cartesian grid and employing boundary layer coupling. It also uses a modal data importer 
and ZAERO 3D spline module to construct structural grids. ZEUS uses central difference with the JST (Jameson- 
Schmidt-Turkel) artificial dissipation scheme for flux construction and Green’s integral boundary layer method for 
turbulence model.”! Aeroelastic analyses in ZEUS require its own input file and a modal solution from a finite element 
solver. In the present work, Nastran is used as a FEM based structural solver for the modal analysis. An aerodynamic 
mesh is generated and the fluid structure interaction is provided in ZEUS to perform the required aeroelastic analyses 
after importing the modal solution. Once the aeroelastic analyses are accomplished for the reference cases, they 
should be repeated as many times as the design samples generated to construct the ROMs. Different flight conditions 
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with various Mach numbers, dynamic pressures, and angle of attack values, for a total of 3 reference cases for flutter 
analyses and 9 reference cases for static aeroelastic analyses are considered. After generating an adequate number 
of sampling for each reference case, the model reduction is performed by using the nonintrusive Polynomial Chaos 
Expansion (PCE) and Proper Orthogonal Decomposition (POD) methods. Both the PCE and POD methods aim to 
ameliorate the computational efficiency. 

Firstly, the static aeroelastic and flutter analyses of the S*T wing were performed for the reference cases by using 
the PCE and POD based ROMs. Then, the ROM results were compared to the initial computational results of the 
high-fidelity models. These ROMs could predict the responses in a small input data range, and in this case both the 
PCE and POD methods give satisfactory results while their accuracies will be thoroughly discussed in the present 
work. 

Secondly, both the PCE and POD methods were used to construct a ROM, which is accurate when large variations 
of input parameters are considered. The flutter behavior of the S*T wing could accurately be predicted over the entire 
flight regime for Mach numbers ranging from 0.60 to 1.20 by the PCE and POD based ROMs. This case will also 
be evaluated in more detail by comparing the flutter results of the ROMs to the initial high-fidelity computational 
analyses. 

Finally, structural and aerodynamic uncertainties, which influence static aeroelastic and flutter computational so- 
lutions, were propagated using ROMs. The uncertainty quantification was performed by using MCS and LHS tech- 
niques, which are frequently used to propagate uncertainties in a probabilistic system. The effects of uncertain input 
variables on the flutter and static aeroelastic responses were investigated. For the static aeroelastic simulations, the un- 
certain input variables examined were angle of attack, elasticity modulus and dynamic pressure. The static aeroelastic 
output variables were aerodynamic lift, drag, moment coefficients and wing tip twist angle. For the flutter simulations, 
the angle of attack and elasticity modulus were the uncertain input variables. The flutter output variables examined 
were flutter frequency, flutter speed and flutter dynamic pressure. 


A. Static Aeroelastic and Flutter Uncertainty Quantification 


In the aeroelastic UQ analysis, we used the distributions for the angle of attack and dynamic pressure as given by 
the AVT-191 task group. Each structural modal solution was used in the aeroelastic analysis and the final results 
were classified according to the given clusters divided for angle of attack and dynamic pressure values. In our study, 
instead of giving a distribution to the modal frequencies, elasticity modulus of Material 21 was distributed by using 
10 samples with a 5% uncertainty around the mean value of 6.39MPa. However, the uncertainty distributions for the 
angle of attack and dynamic pressure were the same as provided by the AVT-191 Task Group. These distributions were 
divided into different clusters and each cluster was defined with a mean value and a standard deviation (see Table 3) 
with a normal distribution assumption. 


Table 3. Clusters for aeroelastic UQ analysis. 


Cluster Angle of attack (degs) Dynamic pressure (psf) 
Number Designation Mean Standard deviation Mean Standard deviation 
1 20 psf Cluster | | 0.204571 0.016328 21.40429 0.182254 
2 20 psf Cluster 2 | 0.826455 0.097638 20.17791 0.171811 
3 35 psf Cluster 0.327333 0.06912 35.55383 0.397882 
4 60 psf Cluster 1 | -0.19118 0.065167 60.9787 1 0.519223 
3 60 psf Cluster 2 | 0.193625 0.043197 60.3455 0.513832 


For the static aeroelastic analysis, the proposed UQ procedure is repeated with three Mach numbers set as M=0.60, 
M=0.85 and M=0.95. The data sampling information in 5 clusters are sequentially run for M=0.60, M=0.85, and 
M=0.95. In Figure 11 to Figure 13, the random aerodynamic lift coefficient C,, drag coefficient Cp, moment coefficient 
Cm and wing tip twist angle values of the S*T wing are plotted against the angle of attack values for the samples in 
each cluster. 

For the flutter analysis, the proposed UQ procedure is similarly repeated with three Mach numbers set as M=0.60, 
M=0.85, and M=0.95. The data samples in the 5 clusters used in our dynamic aeroelastic study are sequentially run for 
M=0.60, M=0.85 and M=0.95. In Figure 14 to Figure 16, the random flutter speed values of the S*T wing are plotted 
against the angle of attack values for the samples in each cluster. 
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Linear aeroelastic analysis methods result in some input-output relationships that are linear and some that are not. 
When the relationship is linear, Gaussian input distributions result in Gaussian output distributions. The nonlinear 
input-output relationships that have been examined here showed that the output data generated due to Gaussian input 
distributions were characterized by beta distributions. The input-output relationships that are not linear, when subject 
to a uniform (constant) input, produced triangular output distributions. Regardless of the output distribution shape, 
maximum likelihood fits to Gaussian distribution parameters were calculated for comparison with other uncertainty 


Figure 11. Cp versus angle of attack for Mach 0.85. 


VI. Discussion 


propagation analyses performed for this project. 


Varying the angle of attack in the linear analysis results in variation only in the trimmed aeroelastic twist angle. 
The static aeroelastic stability and control derivatives are not functions of the angle of attack, and neither are the flutter 
results. The twist angle distributions due to Gaussian-distributed angle of attack variations were demonstrated to be 


directly calculable from the linear relationship properties of the Gaussian statistics. 
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Figure 13. Twist angle versus angle of attack for Mach 0.85. 


Modal analysis showed that variations in the natural frequencies due to the variations in Young’s modulus were 
considerably less than the uncertainty levels that are associated with typical structural dynamic representation. That 
is, the comparison between finite element representation and ground vibration test results is generally greater than the 
observed output variations produced in this study. This indicates that the frequency variations produced by the current 
uncertainty propagation will be considerably lower than should be investigated in a flutter uncertainty study. This 
shortcoming of the current work is expected, as only one region of the FEM was included in the structural property 
input variations. 

On the other hand, the most significant change in the modal frequencies was produced for one of the two modes 
that was shown to be critical for simulating flutter of this configuration. This indicates that the input parameter 
screening process, although performed using static aeroelastic analyses, was effective in choosing a parameter that 
was also significant relative to flutter prediction. In particular, the input parameter screening included evaluating the 
sensitivity of the lift curve slope, Czg, and pitching moment coefficient derivatives due to angle of attack, Cg. While 
these parameters are not direct links to general predictions of flutter, sensitivity of these stability derivatives does 
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Figure 15. Flutter uncertainty quantification for M=0.85. 


provide parameter selection guidance for a configuration where vehicle pitching or wing torsion is involved in the 
flutter mechanism. 

In the second part of the current paper, since the uncertainty quantification using high-fidelity computational mod- 
els requires a very large amount of computational time, ROMs are generated using a high fidelity aeroelastic framework 
based on Euler solver and modal aeroelastic coupling. The static aeroelastic and flutter analyses of the S*T wing are 
performed by using the generated design samples to construct a ROM. The ROM strategy of the present work is based 
on PCE and POD techniques. After the construction of the PCE and POD based ROMs, the reference aeroelastic anal- 
yses are also performed via these models and their accuracies are evaluated. Finally, the structural and aerodynamic 
uncertainties are propagated through the reduced models to investigate the randomnesses in the output parameters. 
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Figure 16. Flutter uncertainty quantification for M=0.95. 


VII. Concluding Remarks 


The influences of the input distribution type and sampling were examined for several cases. It was demonstrated 
that where the input-output relationship is approximately linear, an input distribution shape is reproduced in the output 
distribution shape for both uniform and Gaussian processes. For cases where the relationship is not linear, a Gaussian 
input distribution can produce output parameters that are skewed and better fit with beta distributions. 

Linear analyses are useful for determining which input parameters should have uncertainty distributions assigned. 
Development of statistical distributions to represent those uncertainties relies heavily on observations made in ex- 
perimental data sets. Distributions of the input parameters were determined by examining experimental data and by 
examining the structural model tuning parameters and values that were used by the experiment test team. Parame- 
ter screenings using static aeroelastic stability and control coefficients were performed to determine which material 
property to vary. The flexible-to-rigid ratio of the lift coefficient due to angle of attack and due to the horizontal 
tail produced the same result: the Young’s modulus of the fuselage elements at the wing centerline had the greatest 
influence on these aeroelastic outputs. Uncertainty propagation studies and the parameter screening analyses showed 
that the variation introduced into the lift coefficient due to the horizontal tail is the most sensitive derivative relative to 
Young’s modulus variations. 

The linear analysis method employed was a perturbation analysis and thus can not be used to assess the influence 
of angle of attack uncertainty. This assessment was addressed with higher fidelity modeling, which utilized reduced 
order methods to perform the uncertainty propagation. The results of this study indicate that uncertainty in the Young’s 
modulus associated with the structural elements of the flexible fuselage beam produce significant uncertainties in the 
static aeroelastic and flutter results. 

The results of this study show that the flutter solution is very sensitive to the structural damping. The S*T team 
showed a considerable amount of variation in their attempts to measure the modal damping. These two points lead to 
the undesirable conclusion that the flutter analysis is sensitive to a quantity that is difficult to accurately measure. The 
degree to which the flutter results are sensitive to the structural damping is dependent on, and inversely related to, the 
steepness of the flutter crossing. 
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